R como SIG
Julián Equihua
R como una plataforma SIG
El análisis espacial es una extensión natural a las capacidades de R. Existen múltiples paquetes para manipular y analizar datos espacio-temporales en R. En esta clase vamos a dar una introducción a esto.
A la hora de llevar a cabo visualizaciones y análisis espaciales en R se heredan todos los beneficios que se han venido mencionando en clases sesiones pasadas.
- R permitirá repetir y replicar procesos una ves que se han desarrollado
- Por lo anterior se genera de manera automática investigación transparente
- R es primordialmente un lenguaje de programación para análisis de datos por lo que un sin fin de técnicas estadísticas que no se encuentran tradicionalmente en plataformas SIG están a nuestra disposición
El contra más grande de usar R pos sí solo como SIG
- Las capacidades de visualización de datos espaciales en R, en general, son estáticas. Si se desea explorar un mapa de manera dinámica (Zoom in o Zoom out) se deben de utilizar funcionalidades específicas.
Por esto es recomendable parear R con software SIG especializado, por ejemplo QGIS, que es sumamente poderoso y de código abierto:
Formatos de datos espaciales
En el mundo del análisis espacial existen dos formatos principales:
- El formato vectorial
- El formato raster
El formato vectorial, descripción y ejemplos
El formato vectorial se refiere a estructuras que se conforman por puntos, líneas o polígonos asociados a una tabla de datos. Por ejemplo poner unos ejemplos:
Cada instancia de un crimen asociado a un lugar (un par de coordenadas) se puede asociar a un formato vectorial de puntos. Por ejemplo un mapa de Homicidios dolosos en la CDMX desde Enero del 2013 se podría visualizar de la siguiente manera usando un formato vectorial:
Las carreteras y caminos son candidatos naturales para digitalizarse como líneas.
cada delimitación política de un estado de la república es un polígono, por lo tanto una manera natural de almacenar digitalmente un mapa de los estados es a través de un archivo vectorial. Como se mencionó, estas unidades (polígonos en este caso) pueden estar asociadas a una tabla de datos, por ejemplo, las variables extensión total en Ha y población total.
Un primer objeto espacial en R e introducción a APIs
Casi todas las estructuras de datos vectoriales en R se basan en el paquete. Hoy en día es prácticamente un estándar y la mayoría de los otros paquetes construyen sobre este.
El paquete sp está construído de manera tal que está dividido en funcionalidades que permiten crear y manipular los objetos de la sección anterior: puntos, líneas y polígonos.
Para incentivar el entendimiento de estas estructuras de datos, construyamos un objeto espacial de cero usando el paquete sp.
El geocoding es el proceso computacional de transformar una dirección (texto) a coordenadas geográficas. Existen numerosos sevicios para llevar esto a cabo, el más popular es el de Google Geocoding API:
https://developers.google.com/maps/documentation/geocoding/start
API son las siglas en inglés para Application Programming Interface. Una API no es lo mismo que un servidor remoto, es una parte de un servidor que recibe peticiones y regresa respuestas, esto es; una API es una manera de ofrecer servicios a clientes. En este caso, Google, a través de la API anterior ofrece el servicio de geocoding.
Para usarla se debe generar y activar una key en esta dirección:
https://developers.google.com/maps/documentation/geocoding/get-api-key
Luego conviene almacenarla como un objeto en R.
key_google_geocoding = "AIzaSyCup9GDTLlNdZgc9F_sUZjObp4UynwmNNM"En R existen numerosos paquetes para llevar a cabo consultas de geocoding, a continuación usaremos el paquete “googleway” para obtener datos (entre ellos las coordenadas) de un conjunto de cines cinépolis.
library("googleway")## Warning: package 'googleway' was built under R version 3.4.2resultados_cinepolis = google_places(search_string = "cinepolis,CDMX,Mexico",
place_type="movie_theater",key = key_google_geocoding) # se debe de utilizar la key anterior
names(resultados_cinepolis)## [1] "html_attributions" "next_page_token" "results"
## [4] "status"La función básica “names()” nos regresa los nombres de los elementos de una lista o de las columnas de un data.frame. Usándola sobre el objeto anterior, vimos que en la lista se encuentran los objetos “html_attributions”, “next_page_token”, “results” y “status”.
El objeto “results” es un data.frame que tiene todos los datos sobre las direcciones encontradas al llevar a cabo el proceso de geocoding. De estas vamos a elegir un subconjunto útil.
library("tidyverse")## Loading tidyverse: ggplot2
## Loading tidyverse: tibble
## Loading tidyverse: tidyr
## Loading tidyverse: readr
## Loading tidyverse: purrr
## Loading tidyverse: dplyr## Warning: package 'dplyr' was built under R version 3.4.2## Conflicts with tidy packages ----------------------------------------------## filter(): dplyr, stats
## lag(): dplyr, statscinepolis = data.frame(select(resultados_cinepolis$results,place_id,name,formatted_address,rating),
lat = resultados_cinepolis$results$geometry$location$lat,
lon = resultados_cinepolis$results$geometry$location$lng,
stringsAsFactors = FALSE)El data.frame “cinepolis” todavía no es un objeto espacial. Aún así, utilizando la función “gvisMap()” en el paquete “googleVis” podemos visualizar estos cines de manera dinámica. El parámetro esencial de esta función es el nombre de la columna donde se encuentran las coordenadas de los puntos, en esta columna las coordenadas deben estar en formato texto y con la forma:
lat:lon
construyamos un objeto con esta forma y agreguémoslo como columna a nuestro data.frame cinepolis.
cinepolis$google_latlon = paste0(cinepolis$lat,
rep(":",nrow(cinepolis)),
cinepolis$lon)Ya que este data.frame tiene esta columna en el formato apropiado podemos utilizar la función “gvisMAP()”, el argumento “tipvar” indica qué info estará asociada a los marcadores en el mapa que desplegaremos.
library("googleVis")## Warning: package 'googleVis' was built under R version 3.4.2## Creating a generic function for 'toJSON' from package 'jsonlite' in package 'googleVis'##
## Welcome to googleVis version 0.6.2
##
## Please read Google's Terms of Use
## before you start using the package:
## https://developers.google.com/terms/
##
## Note, the plot method of googleVis will by default use
## the standard browser to display its output.
##
## See the googleVis package vignettes for more details,
## or visit http://github.com/mages/googleVis.
##
## To suppress this message use:
## suppressPackageStartupMessages(library(googleVis))cinepolis_en_google <- gvisMap(cinepolis, locationvar="google_latlon", tipvar="name",
options=list(showTip=TRUE, showLine=F, enableScrollWheel=TRUE,
mapType='satellite', useMapTypeControl=TRUE, width=600,height=400))
# luego las funciones
# plot(cinepolie_en_google) abre una ventana con el mapa
# print(cinepolis_en_google, file="tu_ruta/cinepolis.html") guarda el mapa en un archivo html¿Qué es una proyección cartográfica? - objetos espaciales en R
Ahora usaremos los paquetes “rgdal” y “sp” para convertir el data.frame anterior a uno espacial usando la función “coordinates()” para indicarle a R qué columnas serán sus coordenadas. Esto convierte al objeto “cinepolis” en un data.frame espacial.
library("sp")
library("rgdal")## Warning: package 'rgdal' was built under R version 3.4.2## rgdal: version: 1.2-15, (SVN revision 691)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
## Path to GDAL shared files: D:/software/R/R-3.4.1/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: D:/software/R/R-3.4.1/library/rgdal/proj
## Linking to sp version: 1.2-5coordinates(cinepolis) =~ lon + lat
class(cinepolis)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"Proyectar es llevar a cabo un proceso para hacer que un objeto sea visible sobre otro. En geografía, una proyección se basa en un modelo de la tierra para mapear coordenadas de la misma hacia un plano (x,y).
En este momento nuestra data.frame cinepolis no está proyectada:
summary(cinepolis)## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## lon -99.25832 -99.01790
## lat 19.29530 19.50611
## Is projected: NA
## proj4string : [NA]
## Number of points: 20
## Data attributes:
## place_id name formatted_address rating
## Length:20 Length:20 Length:20 Min. :3.700
## Class :character Class :character Class :character 1st Qu.:4.075
## Mode :character Mode :character Mode :character Median :4.200
## Mean :4.220
## 3rd Qu.:4.400
## Max. :4.600
## google_latlon
## Length:20
## Class :character
## Mode :character
##
##
## Existen muchas proyecciones, la oficial para nuestro país (INEGI) es una basada en un modelo cónico de la tierra. Esta refiere elementos que se hallan en la esfera de la Tierra sobre un cono tangente, usando de vértice el eje que une a los dos polos.
Una proyección particularmente común es WGS84 latitude-longitude, por ejemplo es la utilizada por google maps.
un Coordinate Reference System (CRS) es la combinación de un sistema de coordenadas geográficas y posiblemente de una proyección. Ya que tenemos nuestros coordenadas asignemosle una proyección. Esto se puede hacer de múltiples maneras, las proyecciónes populares tienen asignado un código:
- WGS84 latitude-longitude := CRS(“+init=EPSG:4326”)
O se puede definir parámetro por parámetro:
- WGS84 latitude-longitude := CRS(“+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs”)
Algunos de los parámetros son evidentes, por ejemplo +ellps=WGS84 se refiere a que el elipsoide 3d que se usa para modelar a la tierra es uno llamado WGS84. Sabemos que las coordenadas serán longitud y latitud. Aún así, entender los parámetros de distintas proyecciones escapa el alcance de esta clase.
Lo más importante aquí es saber que si se quieren asociar espacialmente dos conjuntos de información (por ejemplo dos mapas vectoriales), deben encontrarse en la misma proyección.
Asignemos a nuestro data.frame la proyección WGS84 latitude-longitude, esto se puede hacer con la función proj4string()
# asignamos una proyección a estos puntos
proj4string(cinepolis) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")Ahora, este objeto está bien preparado para escribirse en disco y usarse como una verdadera capa SIG.
#writeOGR(cinepolis, "cinepolis.shp", "cinepolis", driver="ESRI Shapefile")Este archivo se puede cargar nuevamente como un data.frame espacial con la función “readOGR()” en el paquete “rgdal”
# ruta_shapefile = "mi_ruta/cinepolis.shp"
# data_frame_espacial = readOGR(ruta_shapefile,ogrListLayers(ruta_shapefile)[1])
Cargar el shapefile generado y con base en él volver a generar el mapa con la función “gVisMap()” nota que esta función recibe un data.frame normal, no uno espacial.
Un data.frame espacial está compuesto por una tabla de datos asociada a una estructura espacial (geometría). Se puede acceder a la tabla de datos haciendo
# data.frame@dataO se puede convertir el data.frame espacial a uno usual con la función “as.data.frame()”.
Por otro lado hay que tener cuidado con los nombres de las columnas puesto en el proceso de escribir un shapefile a disco generalmente se acortan los nombres de columna muy largos.
Cómo usar el DENUE y no morir en el intento
Las funciones de geocoding no están diseñadas para generar datos exhaustivos sobre un negocio en particular. De hecho, el límite de resultados que te puede arrojar una petición en particular es de 20, esto 3 veces utilizando el page_token para un total de 60 resultados (ver la ayuda de la función “google_places()”).
Existen trucos para sobrellevar esta limitante pero en general no es la mejor manera de abordar el problema de ubicar exhaustivamente negocios, por ejemplo haciendo algo como lo que se presenta aquí (aunque usan el lenguaje python se puede implementar en R):
https://iliauk.com/2015/12/18/data-mining-google-places-cafe-nero-example/
En México existe una base de datos decente de negocios. El DENUE, que ahora está asociado a una API y al paquete de R “inegiR” para hacer consultas desarrollados por INEGI mismo.
Para utilizar este paquete de R, análogamente a como se hizo con la API de Google Geocoding se debe generar una key (nota que el DENUE tiene su propio generador de keys):
http://www.beta.inegi.org.mx/servicios/api_biinegi.html
Por supuesto también es buena idea guardarla en un objeto de R.
# asignar key a un objeto
key_denue = "b93b9e6c-5ec5-482d-bce4-03321c402f70"Hay dos maneras de trabajar con el paqueteR para obtener datos del INEGI:
- Una general en la que se hace un parse de un URL provisto por un usuario. Todos los URL’s para cada serie de datos se puede encontrar en el sitio de desarrollo del INEGI:
http://www.beta.inegi.org.mx/servicios/api_biinegi.html
- Una con base en las funciones ya incluídas en el paquete “inegiR”
Por suerte, los datos del DENUE se encuentran en este segundo caso. Por ejemplo la función “denue_inegi()” permite hacer una búsqueda espacial alrededor de un punto de referencia. Por ejemplo alrededor de un cine de nuestra elección.
library("inegiR")## Warning: package 'inegiR' was built under R version 3.4.2# elijamos un cine
cine = as.data.frame(cinepolis[2,])
# luego hacemos la consulta espacial buscando por ejemplo cines cinemex
cinemex_alrededor <- denue_inegi(latitud = cine$lat,
longitud = cine$lon,
token = key_denue,
metros=5000,
keyword="cinemex")
cinemex_alrededor## id Nombre Razon
## 1 6309618 CINEMEX TEZONTLE OPERADORA DE CINEMAS SA DE CV
## 2 6314053 CINEMEX PLAZA ORIENTE OPERADORA DE CINEMAS SA DE CV
## Actividad
## 1 Exhibición de películas y otros materiales audiovisuales
## 2 Exhibición de películas y otros materiales audiovisuales
## Estrato Vialidad Calle NumExterior NumInterior
## 1 31 a 50 personas AVENIDA CANAL DE TEZONTLE 1512
## 2 31 a 50 personas AVENIDA CANAL DE TEZONTLE 1520
## Colonia CP Ubicacion
## 1 DR ALONSO ORTIZ TIRADO 9020 IZTAPALAPA, Iztapalapa, CIUDAD DE MÉXICO
## 2 DR ALFONSO ORTIZ TIRADO 9020 IZTAPALAPA, Iztapalapa, CIUDAD DE MÉXICO
## Tel eMail SitioWeb Tipo Longitud Latitud
## 1 Fijo -99.08240723 19.38199430
## 2 Fijo -99.08169257 19.38182559# podríamos estar interesados también en otros cines menos conocidos
# para estos se puede primero buscar negocios cuyo nombre incluya
# la palabra cine
negocios_alrededor <- denue_inegi(latitud = cine$lat,
longitud = cine$lon,
token = key_denue,
metros=5000,
keyword="cine")
# luego descartamos todos los negocios que no sean cines
cines_alrededor = filter(negocios_alrededor,
Actividad=="Exhibición de películas y otros materiales audiovisuales")
cines_alrededor## id Nombre Razon
## 1 1031205 CINE NACIONAL
## 2 891138 LA CASA DEL CINE MX CONSULADO CINEMATOGRÁFICO, A.C.
## 3 1001689 CINE VENUS DISTRITUR, S.A. DE C.V.
## Actividad Estrato
## 1 Exhibición de películas y otros materiales audiovisuales 6 a 10 personas
## 2 Exhibición de películas y otros materiales audiovisuales 6 a 10 personas
## 3 Exhibición de películas y otros materiales audiovisuales 0 a 5 personas
## Vialidad Calle NumExterior NumInterior Colonia
## 1 AVENIDA FRAY SERVANDO TERESA DE MIER 0 CENTRO
## 2 CALLE REPÚBLICA DE URUGUAY 52 2 CENTRO
## 3 CALLE REPÚBLICA DE CHILE 34 CENTRO
## CP Ubicacion Tel
## 1 06090 CUAUHTÉMOC, Cuauhtémoc, CIUDAD DE MÉXICO 5555224024
## 2 06000 CUAUHTÉMOC, Cuauhtémoc, CIUDAD DE MÉXICO 5555124243
## 3 06020 CUAUHTÉMOC, Cuauhtémoc, CIUDAD DE MÉXICO 5555231099
## eMail SitioWeb Tipo Longitud
## 1 ESTRELLA.MARINLOPEZ@YAHOO.COM.MX Fijo -99.12994101
## 2 CONTACTO@LACASADELCINE.MX WWW.LACASADELCINE.MX Fijo -99.13762812
## 3 ESTRELLA.MARINLOPEZ@YAHOO.COM.MX Fijo -99.13563358
## Latitud
## 1 19.42352399
## 2 19.43067352
## 3 19.43788506Por supuesto lo anterior nos permite automatizar el proceso de búsqueda de otros cines alrededor de cines cinepolis.
Una breve introducción a iteración
Una herramienta útil para evitar escribir más código son los mecanismos de iteración. Estos son útiles cuando se quiere llevar a cabo la misma tarea múltiples veces.
En esta clase aprenderemos sobre un paradigma de iteración: la programación imperativa.
Este paradigma es el más antiguo y permite introducirse fácilmente al tema pues hace que la iteración sea muy explícita. Como desventaja, las estructuras de este paradigma, llamados bucles (en inglés loops), tienden a ser verbosos.
Existe otro paradigma llamado programación funcional que ofrece herramientas para extraer todo el código duplicado para que cada bucle tenga su propia función. Después de un poco de práctica se puede resolver los problemas más comunes en iteración de manera más sencilla, utilizando menos código y por lo mismo cometiendo menos errores. Revisar el paquete “purrr”:
https://cran.r-project.org/web/packages/purrr/index.html
Programación imperativa
Un bucle está compuesto de una secuencia y un cuerpo que se ejecuta con base en esta secuencia.
Una secuencia:
1:5## [1] 1 2 3 4 5Un bucle for puede ejecutar el mismo código (el contenido en su cuerpo) variando un iterador. Por ejemplo el siguiente código nos muestra uno por uno los números en la secuencia misma y en cada una de esas iteraciones nos muestra la palabra “gatito”.
for (i in 1:5) # para la secuancia del 1 al 10
{
print(i)
print("gatito") # ejecutar esto
}## [1] 1
## [1] "gatito"
## [1] 2
## [1] "gatito"
## [1] 3
## [1] "gatito"
## [1] 4
## [1] "gatito"
## [1] 5
## [1] "gatito" Ejercicio: Escribir una bucle que encuentre los cinemex alrededor de cada uno de los cinepolis de nuestro data.frame “cinepolis”, que convierta los data.frames resultantes de esa búsqueda en data.frames espaciales y que los guarde en disco con el nombre del cinepolis respectivo.
Cálculo de distancias entre puntos geográficos
R ofrece paquetes muy convenientes para calcular distancias y tiempos de recorrido entre localidades geográficas. Por ejemplo el paquete “gmapsdistance” permite calcular distancias y tiempos con distintos medios de transporte y partiendo en una fecha:hora determinada (sólo si se brina una key).
Por ejemplo podemos calcular la distancia y tiempo recorrido partiendo ahorita (es el default si no se especifica la fecha:hora de partida) al transportarse entre Toluca y CDMX:
library("gmapsdistance")## Warning: package 'gmapsdistance' was built under R version 3.4.2resultados = gmapsdistance(origin = "Toluca+Estado+De+Mexico+Mexico",
destination = "CDMX+Mexico",
mode = "driving")
resultados## $Time
## [1] 4303
##
## $Distance
## [1] 67749
##
## $Status
## [1] "OK"Esta función también entiende coordenadas lat lon, por ejemplo podemos encontrar la distancia y tiempo de recorrido caminando de la ciudad griega de Marathon a Atenas. Athens.
resultados = gmapsdistance(origin = "38.1621328+24.0029257",
destination = "37.9908372+23.7383394",
mode = "walking")
resultados## $Time
## [1] 30487
##
## $Distance
## [1] 39509
##
## $Status
## [1] "OK"Finalmente, se puede introducir a los argumentos “origin” y “destination” arreglos de localidades y se calcularán las combinaciones pertinentes.
or = c("Washington+DC", "New+York+NY", "Seattle+WA", "Miami+FL")
des = c("Los+Angeles+CA", "Austin+TX", "Chicago+IL", "Philadelphia+PA")
resultados = gmapsdistance(origin = or, destination = des, mode = "bicycling", combinations = "pairwise")
resultados## $Time
## or de Time
## 1 Washington+DC Los+Angeles+CA 858830
## 2 New+York+NY Austin+TX 598038
## 3 Seattle+WA Chicago+IL 675002
## 4 Miami+FL Philadelphia+PA 411413
##
## $Distance
## or de Distance
## 1 Washington+DC Los+Angeles+CA 4582527
## 2 New+York+NY Austin+TX 3138643
## 3 Seattle+WA Chicago+IL 3587447
## 4 Miami+FL Philadelphia+PA 2173558
##
## $Status
## or de status
## 1 Washington+DC Los+Angeles+CA OK
## 2 New+York+NY Austin+TX OK
## 3 Seattle+WA Chicago+IL OK
## 4 Miami+FL Philadelphia+PA OK Ejercicio: Calcular las distancias y tiempo de recorrido desde la CDMX a todos los puertos del país.